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Abstract 

A powerful approach to solve the Coulombic quantum three-body problem is proposed. The 
approach is exponentially convergent and more efficient than the Hyperspherical Coordinate(HC) 
method and the Correlation Function Hyperspherical Harmonic(CFHH) method. This approach is 
numerically competitive with the variational methods, such as that using the Hylleraas-type basis 
functions. Numerical comparisons are made to demonstrate them, by calculating the non-relativistic 
& infinite-nuclear-mass limit of the ground state energy of the helium atom. The exponentially 
convergency of this approach is due to the full matching between the analytical structure of the 
basis functions that I use and the true wave function. This full matching was not reached by almost 
any other methods. For example, the variational method using the Hylleraas-type basis does not 
reflects the logarithmic singularity of the true wave function at the origin as predicted by Bartlett 
and Fock. Two important approaches are proposed in this work to reach this full matching: the 
coordinate transformation method and the asymptotic series method. Besides these, this work 
makes use of the least square method to substitute complicated numerical integrations in solving 
the Schrodinger equation, without much loss of accuracy; this method is routinely used by people to 
fit a theoretical curve with discrete experimental data, but I use it here to simplify the computation. 

PACS number(s): 



1 INTRODUCTION 

Most approximate methods to solve a linear partial differential equation, such as the stationary state 
Schrodinger equation, are actually to choose an A^-dimensional subspace of the infinite-dimensional 
Hilbert space and then to reduce the partial differential equation to N linear algebraic equations 
defined in this subspace. The efficiency of this kind of methods is mainly determined by whether one 
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can use sufficient small N to reach sufficient high accuracy, i.e., make the vector most close to the true 
solution in this subspace sufficiently close to the true solution while keeping the dimension A'' not too 
large to handle. 

Most methods to solve the Coulombic quantum three-body problem belong to this class, except for 
some variational methods that make use of some non-linear variational parameters. The differences 
between different methods of this kind mainly lie in different choices of the subspaces of the Hilbert 
space, i.e., different choices of the basis functions to expand the wave function. 

Theoretically, any discrete and complete set of basis functions may be used to expand the wave 
function, and the convergency is easy to fulfilled. But actually, the convergency is often slow and 
makes sufficient accuracy difficult to achieve. The naive hyperspherical harmonic function method [1- 
3] in solving the Coulombic quantum three-body problem is such an example-this slow convergency 
can be illustrated by an analogous and simple example: to expand the function f{x) = \/l — x"^ 
(— 1 < X < +1) as a series of the Legendre polynomials of x. This series is convergent like N~^, where 
s is a positive constant not large and N is the number of Legendre polynomials involved. The reason 
for this slow convergency is that f{x) is singular at x = ±1 but the Legendre polynomials of x are not. 
I call this the mismatching between the analytical structures of the basis functions (the polynomials 
of x) and f{x). 

The correlation function hyperspherical harmonic(CFHH) method[4] were proposed to overcome 
this difficulty. The spirit of this method can be simply illustrated, still using the above example: 
to divide f{x) by an appropriately selected function (called the correlation function) to cancel the 
low order singularities of f{x) at x = ±1, then to expand the remaining function by the Legendre 
polynomials of x. This time, the series is still convergent as N~^, but s is increased by an amount 
depending on how many orders' singularities have been canceled. 

From this simple discussion one can sec that the singularities of the function f{x) are not completely 
canceled by the correlation function, although more sophisticated correlation function can cancel more 
orders' singularities. 

A very simple approach to totally eliminate the singularity is to make an appropriate coordinate 
transformation, and in the same time thoroughly give up the original hyperspherical harmonic function 
method, not just repair it. For example, for /(x) = Vl — x^, one may write x = sm9, where 
— 7r/2 < < tt/2, then /(x) = cos^ and one can expand /(x) as the series about the Legendre 
polynomials of (2/tt)6. This time the series is factorially convergent. The reason is that the analytical 
structures of /(x) and P;((2/7r)0) match-they are both analytical functions on the whole complex 
plane of 6. 

Another useful approach to solve this problem is to use the asymptotic series. Still considering the 



example f{x) = \/\ — x^, one may write the Taylor series 

=f0 + fix + f2X^ + f3X^ + -- - . 

Of course, this series is slowly convergent near x = ±1. But one can use the following asymptotic 
series to calculate /„ when n is large: 

/„ = ((-1)- + l)(c3/2n-3/' + C5/2n-^/' + cy^n-'/^ + ■■■), 



or, equivalently. 



1/2+L , 
Si 



u = ((-ir + 1) E fs u M ' 



where s = ^, |, |, • • • , \ + L, and s\ = r(s+l). For a given n ^ 1, the error of this formula is minimized 
when L/n ~ 2/3, and the minimum error is about \/-^n~^3~"', which exponentially decreases with 
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n increasing. Using such kind of asymptotic formulae to calculate the high order coefficients of the 
Taylor series, one can expand the singular function f{x) at high precision, with only finite linear 
parameters, /o, • • • , /n and /1/2, • • • , /1/2+L • 

Now I introduce an alternative approach to reduce a differential equation to a given finite di- 
mensional subspace Cj^f of tbe Hilbert space. Here N is the dimension of the subspace. The central 
problem is how to reduce an operator O in the Hilbert space, e.g., the kinetic energy operator or 
the potential energy operator, to an x matrix in the given subspace. For a state ^ £ Cj^f, the 
state = usually ^ jC,j\/. To reduce O into an N x N matrix means to find a state E G jC,_\f to 
approximate '^o- The usual approach to select S is to minimize 

(H-*o,H-«'o) , 

where (, ) is the innerproduct of the Hilbert space. This approach will reduce O to a matrix with 
elements 

Oij = {(pi,0(pj) , 

where (pi G C^f is a set of orthonormal basis in Cj^j-, satisfying (f)j) = Sij, 1 < i,j < N. In numerical 
calculation, the innerproduct is usually computed by numerical integration, which needs sufficient 
accuracy and might be complicated. An alternative approach that does not need these integrations 
is to write the states as wavefunctions under a particular representation(e.g., the space-coordinate 
representation), and then select H to minimize 

J2HXa)-^0{Xa)f , 
a 

where Xa is some sample points in the defining area of the wavefunctions. In order to ensure S to be a 
good approximation of '^o, the sample points should be appropriately chosen. Usually the number of 
the sample points is greater than and approximately proportional to N, and the separation between 



two neighboring sample points should be less than the least quasi-semiwavelength of a wavefunction 
in L]^. 

This alternative approach (I call it the least square method) leads to a reduction of the operator 

O: 

= {4>i,o4>jy , 

where (,)' is a pseudo-innerproduct defined as {(j),ipy = J2a4'*{xa)'>p{xa) for arbitrary cp and ip, and 

01 is a set of pseudo-orthonormal basis in £j\/ satisfying = dij. We find that this approach 
is very similar to the usual one, except that a discrete sum over sample points takes the place of the 
usual innerproduct integration. And there is a great degree of freedom in the selection of the sample 
points. In fact, as soon as the sample points are selected according to the spirit mentioned above, the 
accuracy of the solution of the differential equation usually will not decrease significantly. The major 
factor that determines the accuracy of the solution is the choice of the subspace Cj\/, which has been 
discussed to some extent in previous pages. 

In this work, solving the simplicst quantum three-body problem, the three methods discussed 
above are all used: the coordinate transformation method, the asymptotic series method, and the 
least square method. A high precision is reached for the ground state energy of the ideal helium atom, 
and the solution has also some merit in comparison with the Hyleraas-type variational solution [5, 6]. 
In section 2 the Bartlett-Fock expansion[7,8,9] is studied, in order to reflect the analytical structure of 
the wavefunction near the origin. In this study, the asymptotic series are used to represent the hyper- 
angular dependence of the wavefunction. In section 3 the {u, w) coordinate system is used to study the 
hyper-angular dependence of the wavefunction. This coordinate system cancels the singularity of the 
hyper-angular functions totally. The relationship between this coordinate system and the Hyleraas- 
type variational method is also discussed. The least square method is used to reduce the hyper-angular 
parts of the kinetic energy operator and the potential energy operator to finite-dimensional matrices. 
In section 4 the connection of the outer region solution and the inner region Bartlett-Fock expansion 
is studied, using the least square method. In section 5 the numerical result is presented and compared 
with those of other methods. Some explanations are made. In section 6 some discussions are presented 
and some future developments are pointed out. 

2 BARTLETT-FOCK EXPANSION 

Considering an S state of an ideal helium atom, that is, assuming an infinite massive nucleus and 
infinite light speed, one may write the Schrodinger equation 

-2t{dl + dl + dl + -d,)i; + vij = E^, (1) 



where x = — r^, y = 2rir2cos0i2, z = 2rir2sin^i2, and t = + f"! ~ ^/x^ + + z^. ri 
and r2 are the distances of the electrons from the nucleus, and O12 is the angle formed by the two 
electronic position vectors measured from the nucleus. In this equation, an S state is assumed, so the 
wavefunction ^ is only dependent on ri, r2 and 9 12, or, equivalently, x, y, and z. The atomic unit, 
i.e., h = me = e^/(47r£o) = 1, is assumed throughout this paper. The potential energy is 

V = -l-l^±. (2, 

ri r2 ri2 

where ri2 is the distance between the two electrons. 



''I = V ' = \ ^-7r- ' = \/t-y ■ (3) 



The Bartlett-Fock expansion is 

V' = EV'n,.i'^^' (4) 
n,k 

where n = 0, 1/2, 1, 3/2, 2, • • •, and k = 0, 1, 2, • • •. V'n.fe only depends on the two hyper-angles, say, 
a = x/t and (3 = y/t, and does not depend on the hyper-radius, p = \/t. When k > n, 'ipn,k = 0. 
Using the coordinates t, a, and /3, one may rewrite the Schrodinger equation (1) as 

(9| + ^dt + ^Lo)V = ivt-'/' + Pt'')i^ , (5) 

where p = —E/2, and 

Lo = (1 - a2)a^ - 2a/?aa5;3 + (1 - I3^)d} - Sad^ - Sf^dp ; (6) 

^ 1/2 

^ = , , H — , . (7) 

VlT^ ^ ' 

Substituting eq.(4) into eq.(5), and comparing the corresponding coefficients before t"(lnf)'^, one 
will obtain 

^nV'n,fc + (2n + 2)V'„,fc+l + V'n,fc+2 = Vil^n-^,k + P'^n-l,k , (8) 

where L„ = n(n + 2) + Lq. 

The functions -^^^fc are solved out in the order with n increasing; and for each n, with k decreasing. 
The physical area of (a,/?) is the unit circle: + < 1. And the function ■0„^fe(a,/?) may has 
singularities at a = ±1 and at /? = 1. The singularities are of these kinds: (1 — a)*, (1 + aY, and 
(1 — py, with s = ^,|,|,---. So one may write the Taylor series in the (a, /3) unit circle: 

oo 

V'n,fc(a,/5) = E ^n,k,a,ba''P'' ■ (9) 

a,b=0 

The singularities make the usual cutoff, a + b < Lf + Lg, inappropriate, because the error decreases 
slowly when Lf + Lg increases. But since we have known the forms of the singularities, we can 



write the asymptotic formulae to calculate those high order Taylor coefficients that have important 
contributions: 



r -1 

^ ' S 



^n,k,a,b = E V5n,fe,(,,. [1 + (-1)1 ; (10 - 1) 



. 1 



Eq.(lO-l) is appropriate when a ^ b and a » 1, while eq.(10-2) is appropriate when b ^ a and 6 ^> 1. 
(^) = (s!)/[a!(s — a)!], and s! = r{s + 1). Here I have assumed the state is a spin-singlet, and thus 
V'n,fe(— ")/?) = V'n,fe(")/5)- For a spin-triplet, the factor [1 -|- (—1)"] in eq.(lO-l) should be substituted 
by [1 -(-!)«]. 

In my actual calculation, the (a, b) plane is divided into four areas: 
the finite area: < a,b < Lf and a + b < Lf + Lg (L/ » Lg S> 1), 
the a-asymptotic area: a > Lf and b < Ls, 
the 6- asymptotic area: b > Lf and a < Lg, 
and the cutoff area: the remain area. 

Eq.(lO-l) is used in the a-asymptotic area, and eq.(10-2) is used in the 6-asymptotic area, while 
the contribution from the cutoff area is neglected for it is extremely tiny when L/ S> S> 1. 

In a word, a relevant hyper-angular function is described by a finite set of parameters up to a 
high precision. These parameters are some Taylor coefficients and some asymptotic coefficients. To 
operate with some functions of this kind means to operate with the corresponding sets of parameters. 
The relevant operations arc: addition of two functions-adding the corresponding parameters of the 
two sets; multiplying a function by a constant-multiplying each parameter in the set by the constant; 
multiplying a function by v{oi, /3)(eq.(7))- an appropriate linear transformation of the set of parameters 
of the multiplied function; solving an equation L„/ = g with g known and f unknown-solving a set of 
linear equations about the parameters corresponding to /. Here, I write the relevant linear equations 
corresponding to the equation L„/ = g: 

[n(n + 2) - (a + 6)(a + 6 + 2)]/„,6 + (a + l)(a + 2)/„+2,6 + {b+ 1)(6 + 2)fa,b+2 = 9a,b ; (H - 0) 

[n{n + 2) - (6 + s){b + s + 2)]/b,, + {s + l){2s + 26 + 3)A,s+i + (6 + 1)(6 + 2)fb+2,s = 9b,s ; (U - 1) 
[n{n + 2)-{a + s){a + s + 2)]/„,, + {s + l)(2a + 2s + 3)/„,,+i + (a + l)(a + 2)/„+2,s = ks ■ " 2) 
The detailed order to solve ipn,k is: 

Case 1: n = [n] + ^, where [n] is an integer. In this case, solve V'n.M from eq.(8„ [„]); and then 
solve tpUjln]-! from eq.(8„j„]_i); • • •; at last solve ipn,o from eq.(8„_o)- For each ipn^k^ the order is: first 
solve the asymptotic coefficients, from s = ^tos = Lj — ^; then solve the Taylor coefficients, from 
a + b = Lf + Ls to a + b = 0(i.e.,a = 6 = 0). 



Case 2: n is an integer. In this case, the order is more compUcated, because the operator L„ has 
zero eigenvalue(s) in this case. The order is as following: 

Step 1: set the asymptotic coefHcients and the a + b > n Taylor coefficients of 'tpn,n to zero; 

step 2: k <— n — 1; 

step 3: if < 0, goto step 8; 

step 4: solve the asymptotic coefficients and a + b> n Taylor coefficients oiipn^k^ from eq.(8„_fe), 
in the order analogous to that of case 1. 

step 5: solve the a + b = n Taylor coefficients of ilJn,k+i, from cq.(8n,k)- 

step 6: solve the a + b < n Taylor coefficients of ipn,k+ij from eq.(8„_fe+i), with a + b decreas- 

ing(analogous to case 1) to 0. 
step 7: k k — 1, and goto step 3; 

step 8: set the a + b = n Taylor coefficients of ipn,o with some free parameters; 

step 9: solve the a + b < n Taylor coefficients of V'n.o, from eq.(8n,o), with a + b decreas- 
ing(analogous to case 1) to 0. 

The free parameters in solving eq.(8)(scc step 8 of case 2) are finally determined by the boundary 
condition: ^ ^ 0, when t +oo. In principle, we can use the Bartlett-Fock expansion (eq.(4)) for 
arbitrary t, because it is always convergent. But actually, when t is large, the convcrgcncy is slow 
and there is canceling of large numbers before this convergency is reached, both of which make the 
Bartlett-Fock expansion impractical. So I only use this expansion when t is relatively small(see ref.[15] 
for similarity): s/t < po. 

In atual calculation, I chose Lf = 100, Lg = 20, Lj = 6, rimax = 7.5 (the largest n value of the 
terms in eq.(4) that are not neglected), and po = 0.4, and found that the numerical error for the 
calculation of the inner region < po) wavefunction is no more than a few parts in 10^°. I use 
this method to test the accuracy of the calculation: set E in eq.(8) (note that p = —E/2) equal to 
an initial value (for example, set EinUiai = —2.9037, or set EinUiai = —2.903724377), and use the 
approximate wavefunction ipapp thus obtained to calculate the value {Hil;app)/'4'app, where H is the 
exact Hamiltonian operator, and I find it to be almost equal to the initial value EinUiah with a relative 
error no more than a few parts in 10^^. 

When t is larger, another approach is used: 



3 THE HYPER- ANGULAR DEPENDENCE OF THE WAVEFUNC- 
TION 



We have seen that the hyper-angular dependence of the wavefunction, described as a function of 



(a,/?) for each fixed p = \ft = \jr\ + rf, has singularities at a = ±1 and at /3 = 1. Physically, 
this corresponds to the case that the distance between two of the three particles equals zero. It 
can be proved that, for a spin-singlet, the following coordinate transformation will eliminate these 
singularities thoroughly: 



u 



1 + Q. 1 — a 



Equivalently, 



+ d^—-l,w = Vl^- (12) 



u = 1 , w = — . (13) 

P P 



If the energy-eigenstate ip is symmetric under the exchange of ri and r2(spin-singlct), I believe that, 
for each fixed p, ip is a entire function of (u, w). If the energy-eigenstate ^p is antisymmetric under the 
interchange of ri and r2 (spin-triplet), I believe that, for each fixed p, ip = where ^ is a entire 

function of (n, w) . 

This beautiful characteristic makes it especially appropriate to approximate for each fixed p, 
by an n-order polynomial of {u,w), not by an n-order polynomial of (a,/?). The former expansion, 
a polynomial of {u,w), matches the analytical structure of tp; while the latter one, a polynomial of 
(a, /3), does not. The hyper-spherical harmonic function method belongs to the latter expansion, a 
polynomial of (a, (3). So the hyper-spherical harmonic function expansion does not correctly reflect the 
analytical structure of ^. The slow convergency of the hyper-spherical harmonic function expansion 
is only a consequence of this analytical structure mismatching. 

We expect that the {u, w) polynomial expansion converges factorially to the true wavefunction. It 
is worthful to demonstrate a similar example to illustrate this. Consider a function f{x) = exp(— x), — 
1 < X < -1-1; expand f{x) by Legendre polynomials: f{x) = J2i=o flPl{x)'i it can be proved that the 
error of this formula is of the order l/(2"n!), which factorially approach zero as n increases. 

Using the (p, u, w) coordinates, one can write the Schrodinger equation as: 

-kdl + -d, + ^)i^ + -i^ = Ei., (14) 
2 ^ p p 

where Lq and C are the hyper-angular parts of the kinetic energy and the potential energy, respectively. 



u{2 + u) w u{2 + u) w 

(15) 

C = -lll±4 + i. (16) 

u{2 + u) w ^ ^ 



The physical area V of (u, w) is: 
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In this figure, point A corresponds to the coincidence of the two electrons, and point B corresponds 
to the coincidence of the nucleus and one electron. 

For a spin-singlet, we can use an n-order polynomial of (n, w) to approximate V'- The coefficients 
of this polynomial are functions of p. Denote by Cj\f the set of all the polynomials of (u, w) with order 
no more than n. Here, N = [n + l){n + 2)/2 is the dimension. In the physical area P, I choose a set 
of points as sample points: 



/7t ("2 + 0-5) 

i„ = (^/2-l)-[(^/2-l)-mK)]^^i+^ , (18) 



Wa = ^ , (17) 

n2 



where m{w) is the minimum physical u value for a w value. m{w) = \Jl — — 1, if w < 1; and 
m{w) = w — 1, if w > 1. a = (oi, 02), and < ai < ni, < 02 < n2. I chose ni = 77-2 = 2n, so there 
are altogether 4n^ sample points. These sample points define a pseudo-innerproduct. I constructed 
a set of pseudo-orthonormal basis in Cj\f, by using the Schmidt orthogonalization method, and then 
reduce the operators Lq and C to x matrices under this basis, using the method introduced in 
section 1. 

4 CONNECTION OF THE INNER REGION AND THE OUTER 
REGION WAVEFUNCTIONS 

In the area p < po(iiiiier region), the Bartlett-Fock expansion is used. In the area p > po/2(outer 
region), tp is approximated by a vector in Cj\f for each given p, and the partial derivatives with 
respect to p are substituted by optimized variable-order and variable-step differences, which requires 
the selection of a discrete set of p values. The overlap region of the inner region and the outer region 
ensures the natural connection of the derivative of ip, as well as the connection of ip itself. The 
connection is performed by using the least square method: for a polynomial of {u,w) at p = po, 
appropriately choose the values of the free parameters of the solution of eq.(8) (see section 2) so 



that the sum of the squares of the differences of the the inner region solution and the outer region 
polynomial at the sample points is minimized. This defines a linear transformation to calculate the 
values of those free parameters from the given polynomial. When the values of these free parameters 
are determined, one can calculate the values of ip in the region pq/2 < p < pQ, using the Bartlett-Fock 
solution, and further use these ip values to construct polynomials of('u, w) at po/2 < p < Po (according 
to the law of least square), and then use these polynomials in the difference calculation of the partial 
derivative of V' with respect to p at p > pq. At a sufficient large value p = pi, the first-class boundary 
condition is exerted; of course, future development may substitute this by a connection with the long 
range asymptotic solution of V'. 

At last, the whole Schrodinger equation is reduced to an eigen-problem of a finite-dimensional 
matrix. The dimension of the matrix is Np x N, where Np is the number of free p nodes used in 
discretizing the partial derivatives with respect to p, and N is the number of independent hyper- 
angular polynomials used. Note that the energy value should be used in solving eq.(8), but it is 
unknown. The actual calculation is thus an iteration process: choose an initial value of Eq to solve 
eq.(8) and form the Np x dimensional matrix, and calculate the eigenvalue of this matrix to get a 
new value Ei, etc.. The final result is the fixed point of this iteration process. In actual calculation, 
I found that the convergency of this iteration process is very rapid if po is relatively small. Choosing 
Po = 0.4, I found that each step of iteration cause the difference between the eigenvalue of the matrix 
and the fixed point decrease by about (—160) times, when calculating the ground state. 

5 NUMERICAL RESULT AND COMPARISONS 

Using 20 independent Bartlett-Fock series(up to the t^-^ term in eq.(4), neglecting higher order terms), 
choosing n = 10 (so that N = 66), choosing Np = 40, with po = 0-4 and pi = 11.32, and with the 
discrete values of p equal to 0.4/1.2^ 0.4/1.2^, 0.4/1.2, 0.4, 0.4x 1.2, 0.4x 1.2^, 0.4x 1.2^ • • • , 0.4x 1.2^ = 
1.7199, 0.4 x 1.2*^ + 0.3, 0.4 x 1.2^ + 0.6, 0.4 x 1.2^ + 0.9, • • • , 0.4 x 1.2^ + 9.3, and 0.4 x 1.2*^+9.6 = 11.32 
(the first three points are for the natural connection of the derivative of tp, the last point is for the 
first-class boundary condition, and the remained 40 points are free nodes), and discretizing the partial 
derivatives with respect to p according to the complex-plane-division rule(that is: when calculating 
the partial derivatives with respect to p at p = Z, use and only use those node points satisfying p> 1/2 
in the difference format, because the point p = is the singular point), I obtained the result for the 
ground state energy of the ideal helium atom: 

E = -2.9037243738 , (19) 



compared with the accurate value: 

E = -2.9037243770 



(20) 



So the relative error of the result (19) is about 1.1 x 10^'^. Since my method is not a variational method, 
the error of the approximate wavefunction that I obtained should be of a similar order of magnitude, 
so if one calculate the expectation value of the Hamiltonian under this approximate wavefunction, the 
accuracy of the energy will be further raised by several orders of magnitude. 

The result (19) is much more accurate than the result of ref. [10]:— 2.90359, which used the hyper- 
spherical coordinate method. In ref. [10], the quantum numbers (Zl, Z2) (angular momenta of the two 
electrons) are used and a cutoff for them is made; this cutoff does not correctly reflect the analytical 
structure of tp at ru = (equivalently /3 = 1). This is the major reason causing the inaccuracy of the 
result of ref. [10]. 

It is also worthful to compare my result with that of ref. [4], in which the correlation function 
hyper-spherical harmonic method is used. Note that the result (19) is obtained by using a set of 
N = 66 hyper-radius-dependent coefficients to expand the wavefunction. For a similar size in ref. [4], 
N=64, the result is -2.903724300, with relative error about 26.5 x lO"'^. When N=169, the result 
of ref. [4] is —2.903724368, with relative error about 3.1 x 10~^. Apparently my method converges 
more rapidly than that of ref. [4]. The major reason is that the correlation function hyper-spherical 
harmonic method does not cancel the singularities totally — there is still some discontinuity for the 
higher order derivatives, although the low order singularities, which trouble the naive hyper spherical 
harmonic method, are canceled by the correlation function. 

6 CONCLUSIONS, DISCUSSIONS AND FUTURE DEVELOPMENTS 

In conclusion, there are several important ideas in my work that should be emphasized: first, I use 
the asymptotic series to compute the Bartlett-Fock series up to a high precision, with error no more 
than, for example, a few parts in 10^°. Second, I propose an alternative coordinate system, the {u,w) 
system, in which the hyper-angular singularities are thoroughly eliminated, which renders a factorial 
convergency for the expansion of the hyper-angular function. Third, I make use of the least square 
method to reduce an operator(infinite dimensional matrix) to a finite dimensional matrix in a finite 
dimensional subspace of the Hilbert space and to connect the solutions in different regions, avoiding 
complicated numerical integrations, without much loss of the accuracy for the solution. Fourth, 
the optimized difference format — the complex plane division rule — is used to discretize the partial 
derivatives of the wavefunction with respect to p. I calculated the ground state energy of an ideal 
helium atom concretely and obtained a very high precision, demonstrating that my method is superior 
to many other methods and competitive with any sophisticated methods. 

About the analytical structure of the stationary wavefunction: 1. there are logarithmic singularities 
at p = 0, in the forms of p'^{h\p)^; 2. for a given if) (for a spin-singlet) or ^^/[{ri — r2)/p](for a 
spin-triplet) has no singularity, as a function of {u,w). 



Here, I must mention the well known variational method based on the Hyleraas-type functions, 
because it also satisfies the second characteristic of the wavefunction mentioned in the above paragraph. 
One can see this by a simple derivation. The Hyleraas-type function is a entire function of ri, r2 and 
ri2, or equivalently, a entire function of ri + r2, ri — r2, and ri2- For a fixed p, one can substitute 
(ri — r2)^ in this function by — (ri + r2)^, so that, for fixed p, the function is a entire function 
of ri + r2 and ri2 for spin-singlet, or such kind of entire function times a common factor ri — r2 for 
spin-triplet. Equivalently, for fixed p, the Hyleraas-type function is a entire function of (it, (spin- 
singlet) or such kind of entire function times ^''-^^"^ (spin-triplet). This characteristic is one of the most 
important reasons that account for the high accuracy of the Hyleraas-type variational method. 

But this variational method also has its shortcoming: the Hyleraas-type function does not reflect 
the logarithmic singularities with respect to p. So, although this method has high precision for the 
energy levels, the approximate wavefunctions that it renders may deviate significantly from the true 
wavefunctions near the origin. See ref.[4,13,14] for detailed discussions. 

A central idea of this paper is: devising the calculation method according to the analytical structure 
of the true solution. The [u, w) coordinates, the Bartlett-Fock expansion and the asymptotic series 
approach to compute this expansion, and the complex-plane-division rule in calculating the partial 
derivatives with respect to p, all reflect this central idea. The basic principle that ensures high 
numerical precision is just this idea. 

This preliminary work is incomplete in the following aspects: 

First, how to prove that ^(for spin-singlet, or tp/[{ri — r2)/p] for spin-triplet) has no singularity for 
fixed p, as a function of {u,w)'? Note that if this function still has singularities outside of the physical 
area P(see previous figure), the convergency of the expansion of the hyper-angular function will be 
only exponential, not factorial. Of course, even if such kind of singularities do exist, my method will 
still converge more rapidly than the correlation function hyperspherical harmonic method, because 
the latter method only converges like N~p, slower than exp(— 7\/iV). The rapid convergency of my 
method make me guess that such kind of singularities do not exist. 

Second, the asymptotic behavior of the wavefunction, when one electron is far away from the 
nucleus, is not studied in this work. This problem will be important when the highly excited states 
and the scattering states are studied, a topic that will become my next object. 

Third, how to use the ideas proposed in this work to study a helium atom with finite nuclear mass? 
Besides this, the relativistic and QED corrections must be calculated, if one want to obtain a result 
comparable with high-precision experiments. 

Fourth, I have focused on the S states till now. When the total angular momentum is not zero, 
there might be more than one distance-dependent functions (see, for example, ref.[ll]). I believe 
that some important analytical structures of the S states studied in this work are also valid for those 
functions. 



Surely, some important aspects of this work will also play an important role in the highly excited 
states and the scattering states: the logarithmic singularities about p and the method to compute the 
Bartlett-Fock expansion, the non-singularity with respect to the coordinates (u, w) , and the technique 
to connect solutions of different regions, etc.. They can be applied to the study of the highly excited 
states and the scattering states. 
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